options(stringsAsFactors=F)
library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(tidyr)
library(ggrastr)

##########################################################################################

option_list <- list(
  make_option(c("--all_data_by_genes_IGC"),type = "character"),
  make_option(c("--all_data_by_genes_DGC"),type = "character"),
  make_option(c("--amp_gene"),type = "character"),
  make_option(c("--del_gene"),type = "character"),
  make_option(c("--express_file"),type = "character"),
  make_option(c("--dat_info_multi"),type = "character"),
  make_option(c("--dat_info_TCGA"),type = "character"),
  make_option(c("--out_path"), type = "character"), 
  make_option(c("--Class"), type = "character")
)

if(1!=1){
  
"all_data_by_genes.txt"
"CIN_IGCcompareDGC_Driver.copy_number_del_merge.tsv"  
"CombineTMM.DNAUse.NJMU_TCGA.MergeMutiSample.tsv"
"Combine_CIN_sampleInfo.tsv"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

cnv_file1 <- opt$all_data_by_genes_IGC
cnv_file2 <- opt$all_data_by_genes_DGC
amp_gene <- opt$amp_gene
del_gene <- opt$del_gene
express_file <- opt$express_file
dat_info_multi <- opt$dat_info_multi
dat_info_TCGA <- opt$dat_info_TCGA
out_path <- opt$out_path
Class <- opt$Class

###########################################################################################
##画图函数
plotCor <- function(gene = gene , dat_use = dat_use , dat_expression = dat_expression , image_path = image_path,i=i,Class=Class,type=type){
  
  dir.create(image_path , recursive = T)
  
  
  
  title <- paste0(
    "Gene : " , gene 
  )
  
  ## TMM
  tmp_exp <- subset( dat_expression , rownames(dat_expression) == gene)
  dat_use_tmp <- subset( dat_use , rownames(dat_use) == gene)
  ##倒置dat_use和tmp_use
  dat_use_tmp <- data.frame(t(dat_use_tmp))
  dat_use_tmp$sample <- rownames(dat_use_tmp)
  colnames(dat_use_tmp) <- c("copynumber","sample")
  tmp_exp <- data.frame(t(tmp_exp))
  tmp_exp$sample <- rownames(tmp_exp)
  colnames(tmp_exp) <- c("tpm","sample")
  tmp_dat_cor <- merge(dat_use_tmp,tmp_exp,by="sample")
  tmp_dat_cor <- tmp_dat_cor[,-1]
  
  ## 该样本的表达和干细胞评分的关系
  ## TPM
  plot <- ggplot( tmp_dat_cor , aes( x = copynumber , y = tpm ) ) +
    geom_point_rast() + 
    geom_smooth(method = 'lm') +
    ylab("log2(TMM+1)")+
    xlab(paste0("log2(copy number/2) of ",gene))+
    stat_cor(data=tmp_dat_cor, method = "pearson") +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
          legend.position ='none',
          legend.title = element_blank() ,
          panel.grid.major=element_line(colour=NA),
          plot.title = element_text(size = 8,color="black",face='bold'),
          legend.text = element_text(size = 10,color="black",face='bold'),
          axis.text.y = element_text(size = 10,color="black",face='bold'),
          axis.title.x = element_text(size = 10,color="black",face='bold'),
          axis.title.y = element_text(size = 10,color="black",face='bold'),
          axis.ticks.x = element_blank(),
          axis.text.x = element_text(size = 10,color="black",face='bold') ,
          axis.line = element_line(size = 0.5)) 
  out_name <- paste0(image_path,"/",i,"_",gene,"_cor_CnvExpressionOf_",Class,"_",type,".pdf")
  ggsave(file=out_name,plot=plot,width=4,height=4)
}

###########################################################################################
## 导入数据
## 需要cnv数据、amp或del基因的列表、表达矩阵数据
cnv_gene1 <- data.frame(fread(cnv_file1))
cnv_gene2 <- data.frame(fread(cnv_file2))
cnv_gene <- merge(cnv_gene1,cnv_gene2[,c(1,4:ncol(cnv_gene2))],by="Gene.Symbol")

amp_gene <- data.frame(fread(amp_gene))
del_gene <- data.frame(fread(del_gene))
express <- data.frame(fread(express_file))

dat_info_multi <- data.frame(fread(dat_info_multi))
dat_info_TCGA <- data.frame(fread(dat_info_TCGA))

dat_info_TCGA <- subset(dat_info_TCGA,dat_info_TCGA$From=="TCGA" & dat_info_TCGA$Molecular.subtype=="CIN")
dat_info_TCGA$use_sample <- dat_info_TCGA$Tumor
dat_info_TCGA$ID <- dat_info_TCGA$Tumor
dat_info_TCGA$use_sample <- gsub("-",".",dat_info_TCGA$use_sample)
dat_info_TCGA$ID <- gsub("-",".",dat_info_TCGA$use_sample)
dat_info_TCGA$Normal <- dat_info_TCGA$Tumor

dat_info_multi$use_sample <- paste0(dat_info_multi$Tumor,"_",dat_info_multi$Normal)
dat_info_multi <- dat_info_multi[,c("ID","use_sample","Normal","Class.x")]
colnames(dat_info_multi)[4] <- "Class"
dat_info <- rbind(dat_info_multi[,c("ID","use_sample","Normal","Class")],dat_info_TCGA[,c("ID","use_sample","Normal","Class")])
dat_info$sample <- paste0(dat_info$ID,"_",dat_info$Class)
dat <- dat_info
dat$Class <- "all"
dat_info <- rbind(dat_info,dat)

###########################################################################################
##提取amp的基因，并整理cnv表格
for(i in c("amp","del")){
  if(i=="amp"){
    cnv_gene_use <- subset(cnv_gene,Gene.Symbol %in% unique(amp_gene$Gene_Symbol[amp_gene$p < 0.05]))
    }else if(i=="del"){
    cnv_gene_use <- subset(cnv_gene,Gene.Symbol %in% unique(del_gene$Gene_Symbol[del_gene$p < 0.05]))  
    }
  if(nrow(cnv_gene_use)>0){
    rownames(cnv_gene_use) <- cnv_gene_use$Gene.Symbol
    cnv_gene_use <- cnv_gene_use[,-c(1,2,3)]
    ##整合express和cnv_gene_use
    ##现在转录组和CNV的命名方式不一样，需要转成一致的形式
    corTab_all <- data.frame()
    for(type in c("IGC","DGC","all")){
      print(type)
      dat_info_use <- subset(dat_info,dat_info$Class==type)
      express_use <- express[,which(colnames(express) %in% c("gene_id",dat_info_use$sample))]

      dat_info_express <- subset(dat_info_use,dat_info_use$sample %in% colnames(express_use))

      cnv_gene_use_final <- cnv_gene_use[,colnames(cnv_gene_use) %in% dat_info_express$use_sample]
      colnames(cnv_gene_use_final) <- dat_info_express$sample[match(colnames(cnv_gene_use_final),dat_info_express$use_sample)]

      express_use <- express_use[!duplicated(express_use$gene_id),]
      express_use <- subset(express_use,gene_id %in% rownames(cnv_gene_use_final))
      rownames(express_use) <- express_use[,1]
      if(nrow(express_use)>0){
        use_col <- subset(colnames(express_use),colnames(express_use) %in% colnames(cnv_gene_use_final))
        express_use <- express_use[,colnames(express_use) %in% use_col]
        express_use <- log2(express_use+1)
        cnv_gene_use_final <- subset(cnv_gene_use_final,rownames(cnv_gene_use_final) %in% rownames(express_use))
        cnv_gene_use_final <- cnv_gene_use_final[,colnames(cnv_gene_use_final) %in% colnames(express_use)]
        ##把他们的方向调整为一致的方向
        cnv_gene_use_final <- cnv_gene_use_final[ , order(colnames(cnv_gene_use_final), decreasing = FALSE)]
        express_use <- express_use[ , order(colnames(express_use), decreasing = FALSE)]

        ##计算相关性（cnv与表达之间的相关性）
        cna_subset <- cnv_gene_use_final
        expr_subset <- express_use
        corTab <- NULL
        for (j in rownames(cna_subset) ) {
        tmp1 <- as.numeric(cna_subset[j,]) # 相应的CNA数据
        tmp2 <- as.numeric(expr_subset[j,]) # 相应的表达数据
        cor.res <- cor.test(tmp2,tmp1, method = "pearson") # pearson相关性   #Spearman
        #将计算的相关性矩阵合并起来
        corTab <- rbind.data.frame(corTab,
                                   data.frame(gene = j,
                                              Correlation = ifelse(is.na(cor.res$estimate), 0, cor.res$estimate),
                                              Pvalue = ifelse(is.na(cor.res$p.value), 1, cor.res$p.value),
                                              Class = type,
                                              use_nums =length(tmp1),
                                              stringsAsFactors = F),stringsAsFactors = F)
        }

        #out_name <- paste0(out_path,"/",i,"/Gene/cor/cor_CnvExpression_",Class,"_",type,".tsv")
        #write.table(corTab,out_name,sep="\t",row.names=F,quote=F)
        #corTab_all <- rbind(corTab,corTab_all)
        ##相关性大于0.3并且p值小于0.05的才绘图
        if(length(unique(corTab$gene[corTab$Correlation >0.3 & corTab$Pvalue < 0.05])) > 0 ){
          for( gene in unique(corTab$gene[corTab$Correlation >0.3 & corTab$Pvalue < 0.05])){
            plotCor(gene = gene , dat_use = cnv_gene_use_final , dat_expression = express_use , image_path = out_path,i=i,Class=Class,type=type)
          } 
        }
      }
      out_name <- paste0(out_path,"/",i,"_cor_CnvExpression_",Class,".tsv")
      write.table(corTab_all,out_name,sep="\t",row.names=F,quote=F)
    }
  }
}









